#include "stdafx.h"
#include "hnum_cblas.h"

/*  -- LAPACK auxiliary routine (version 2.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       February 29, 1992   

    Purpose   
    =======   

    DLARAN returns a random real number from a uniform (0,1)   
    distribution.   

    Arguments   
    =========   

    ISEED   (input/output) INT array, dimension (4)   
            On entry, the seed of the random number generator; the array 
  
            elements must be between 0 and 4095, and ISEED(4) must be   
            odd.   
            On exit, the seed is updated.   

    Further Details   
    ===============   

    This routine uses a multiplicative congruential method with modulus   
    2**48 and multiplier 33952834046453 (see G.S.Fishman,   
    'Multiplicative congruential random number generators with modulus   
    2**b: an exhaustive analysis for b = 32 and a partial analysis for   
    b = 48', Math. Comp. 189, pp 331-344, 1990).   

    48-bit integers are stored in 4 integer array elements with 12 bits   
    per element. Hence the routine is portable across machines with   
    integers of 32 bits or more.   

    ===================================================================== 
*/
// ===================================================================== 
// Espen Harlinn:
//   This implementation is a slightly modified version version
//   of the code shipped with the SuperLU library, as it makes 
//   no sense to decrement the iseed pointer, and then add 1 to
//   the array offsets to adjust for the decrement.
// ===================================================================== 

namespace harlinn
{
    namespace numerics
    {
        namespace SuperLU
        {
            double dlaran_(int *iseed)
            {
                // multiply the seed by the multiplier modulo 2**48 
                int it4 = iseed[3] * 2549;
                int it3 = it4 / 4096;
                it4 -= it3 << 12;
                it3 = it3 + iseed[2] * 2549 + iseed[3] * 2508;
                int it2 = it3 / 4096;
                it3 -= it2 << 12;
                it2 = it2 + iseed[1] * 2549 + iseed[2] * 2508 + iseed[3] * 322;
                int it1 = it2 / 4096;
                it2 -= it1 << 12;
                it1 = it1 + iseed[0] * 2549 + iseed[1] * 2508 + iseed[2] * 322 + iseed[3] * 494;
                it1 %= 4096;

                // return updated seed 

                iseed[0] = it1;
                iseed[1] = it2;
                iseed[2] = it3;
                iseed[3] = it4;

                // convert 48-bit integer to a real number in the interval (0,1) 
                return ((double) it1 + ((double) it2 + ((double) it3 + (double) it4 * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4;

            } 

        };
    };
};
